#include "ope_term.h"

int Ope_term::phase(int n) { return n % 2 ? -1 : 1; }

int Ope_term::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}

int Ope_term::Tri_state(const State_Jacobi &state) {
  int unit = 1;
  int l12_2 = 2 * state.l12;
  int l3_2 = 2 * state.l3;
  if (Tri(l12_2, state.S12_2, state.J12_2) != 1) {
    return -1;
  }
  if (Tri(l3_2, unit, state.j3_2) != 1) {
    return -1;
  }
  if (Tri(state.J12_2, state.j3_2, state.J_2) != 1) {
    return -1;
  }
  if (Tri(state.T12_2, unit, state.T_2) != 1) {
    return -1;
  }
  return 1;
}

void Ope_term::init() {
  pf_ope = -1.0 * LECs.LEC_CD * LECs.LEC_gA * pow(hbarc, 6) /
           (pow(F_pi, 4) * Cut_LambdaX * 4 * std::sqrt(3.0) * pow(Pi, 4));
  // pf_ope= pf_ope * sqrt(3) * 61528.9083889; // sqrt(3)*(2*pi)^6
}

double Ope_term::cal(const State_Jacobi &final_state, const State_Jacobi &init_state) {

  double result = 0.0;
  if ((Tri_state(final_state) != 1) || (Tri_state(init_state) != 1)) {
    cout << "can`t contain Tri condition in OPE_term.cal";
    exit(0);
    return 0.0;
  }
  // !!!! should n't have E conservation !!!!
  // energy conservation
  //if (final_state.E != init_state.E)
  //  return result;

  if ((final_state.J_2 != init_state.J_2) ||
      (final_state.T_2 != init_state.T_2))
    return result;
  if ((final_state.l12 != 0) || (init_state.l12 != 0))
    return result;
  if ((final_state.J12_2 != final_state.S12_2) ||
      (init_state.J12_2 != init_state.S12_2))
    return result;
  // parity conservation
  if ((final_state.l3 % 2) != (init_state.l3 % 2))
    return result;

  //cout << "basis good ~~" << endl;
  int phase_t = init_state.T_2 + 3 + final_state.j3_2 + init_state.J_2 +
                3 * final_state.l3 + init_state.l3;
  int phase_ope = phase(phase_t / 2);

  /*int phase_test_t =
      init_state.n12 + init_state.n3 + final_state.n12 + final_state.n3;
  phase_ope = phase_ope * phase(phase_test_t);
  */

  double part1 = (final_state.T12_2 + 1) * (init_state.T12_2 + 1) *
                 (final_state.J12_2 + 1) * (init_state.J12_2 + 1) *
                 (final_state.j3_2 + 1) * (init_state.j3_2 + 1);
  part1 = std::sqrt(1.0 * part1);
  double part2 =
      sixJ.quick(final_state.T12_2, init_state.T12_2, 2, 1, 1, 1) *
      sixJ.quick(final_state.T12_2, init_state.T12_2, 2, 1, 1, init_state.T_2);
  if (part2 == 0.0)
    return result;
  double part3 = sixJ.quick(final_state.J12_2, final_state.j3_2, init_state.J_2,
                            init_state.j3_2, init_state.J12_2, 2) *
                 sixJ.quick(1, 1, 2, init_state.J12_2, final_state.J12_2, 1);
  if (part3 == 0.0)
    return result;
  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  if (part1 * part2 * part3 != 0) {
    //cout<<"good ?"<<endl;
    double t_A = term_A(final_state, init_state);
    //cout<<"term A = "<<t_A<<endl;
    double t_B = term_B2(final_state, init_state);
    //cout<<"term B = "<<t_B<<endl;
    result = pf_ope * phase_ope * part1 * part2 * part3 * (t_A + t_B);
    /*
    cout << "check: pf = " << pf_ope << "\tphase = " << phase_ope
         << "\tpart123 =" << part1 * part2 * part3 << endl;
    cout << " term A" << t_A << endl;
    cout << " term B " << t_B << endl;
    */
  }
  return result;
}

double Ope_term::term_A(const State_Jacobi &final_state, const State_Jacobi &init_state) {

  if (final_state.l3 != init_state.l3)
    return 0.0;
  int phase_t = 1 + final_state.j3_2 + 2 * init_state.l3;
  int phase_A = phase(phase_t / 2);
  double part1 =
      1 / 3.0 *
      sixJ.quick(1, 1, 2, init_state.j3_2, final_state.j3_2, 2 * init_state.l3);
  if (part1 == 0.0)
    return 0.0;
  // double part2 = PFK_ope_cal0(final_state, init_state);

  int num_p = wave.nl_vec[final_state.n12][final_state.n3][0][final_state.l3];

  int num = wave.nl_vec[init_state.n12][init_state.n3][0][init_state.l3];

  double part2 = wave.PFKglbar0_vec[num_p][num];

  double result = phase_A * part1 * part2;

  return result;
}

double Ope_term::term_B(const State_Jacobi &final_state, const State_Jacobi &init_state) {

  int n12 = init_state.n12;
  int n3 = init_state.n3;
  int n12p = final_state.n12;
  int n3p = final_state.n3;

  int l3 = init_state.l3;
  int l3p = final_state.l3;
  int phase_B = phase(l3 + l3p);
  double factor_B = 1.0 / (2.0 * std::sqrt(30));

  double part1 = 20 * ninJ.quick(2 * final_state.l3, 2 * init_state.l3, 4, 1, 1,
                                 2, final_state.j3_2, init_state.j3_2, 2);
  if (part1 == 0.0)
    return 0.0;
  // lambda1 = 0; lambda2 = 2
  double part2_0 = phase_B * factor_B * std::sqrt(2.0 * l3p + 1) *
                   cg.quick(2 * l3p, 4, 2 * l3, 0, 0, 0);
  if (part2_0 != 0.0) {
    // part2_0 = part2_0 * PFK_ope_cal2(final_state,init_state,0,2,l3p);
    int num_p = wave.nl_vec[n12p][n3p][0][l3p];
    int num = wave.nl_vec[n12][n3][0][l3];
    part2_0 = part2_0 * wave.PFKglbar2_vec[num_p][num][0][l3p];
  }

  // lambda1 = 1; lambda2 = 1;
  double part2_1 = 0.0;
  int l_min = std::max(std::abs(l3 - 1), std::abs(l3p - 1));
  int l_max = std::min(l3 + 1, l3p + 1);
  if (l_min <= l_max) {
    double part2_1t = 0.0;
    double part2_1tt = 0.0;
    for (int lbar = l_min; lbar <= l_max; lbar++) {
      part2_1tt = 0.5 * (2 * lbar + 1) *
                  cg.quick(2 * lbar, 2, 2 * l3p, 0, 0, 0) *
                  cg.quick(2 * lbar, 2, 2 * l3, 0, 0, 0) *
                  sixJ.quick(2, 2, 4, 2 * l3, 2 * l3p, 2 * lbar);
      if (part2_1tt != 0) {
        // part2_1tt = part2_1tt *
        // PFK_ope_cal2(final_state,init_state,1,1,lbar);
        int num_p = wave.nl_vec[n12p][n3p][0][l3p];
        int num = wave.nl_vec[n12][n3][0][l3];
        part2_1tt = part2_1tt * wave.PFKglbar2_vec[num_p][num][1][lbar];
      }
      part2_1t += part2_1tt;
    }
    part2_1 = part2_1t;
  }

  // lambda1 = 2; lambda2 = 0
  double part2_2 = phase_B * factor_B * std::sqrt(2.0 * l3 + 1) *
                   cg.quick(2 * l3, 4, 2 * l3p, 0, 0, 0);
  if (part2_2 != 0.0) {
    // part2_2 = part2_2 * PFK_ope_cal2(final_state,init_state,2,0,l3);
    int num_p = wave.nl_vec[n12p][n3p][0][l3p];
    int num = wave.nl_vec[n12][n3][0][l3];
    part2_2 = part2_2 * wave.PFKglbar2_vec[num_p][num][2][l3];
  }

  // cout<<" term B :  part1 = "<<part1<<"  part2_0 = "<<part2_0<<"  part2_1 =
  // "<<part2_1<<"  part2_2 = "<<part2_2<<endl;
  double result = part1 * (part2_0 + part2_1 + part2_2);
  return result;
}

double Ope_term::term_B2(const State_Jacobi &final_state, const State_Jacobi &init_state) {

  int n12 = init_state.n12;
  int n3 = init_state.n3;
  int n12p = final_state.n12;
  int n3p = final_state.n3;

  int l3 = init_state.l3;
  int l3p = final_state.l3;
  int phase_B = phase(l3 + l3p);

  double part1 =
      1.82574185835 * ninJ.quick(2 * l3p, 2 * l3, 4, 1, 1, 2, final_state.j3_2,
                                 init_state.j3_2, 2);
  if (part1 == 0.0)
    return 0.0;
  // lambda1 = 0; lambda2 = 2
  double part2_0 = phase_B * std::sqrt(2.0 * l3p + 1) *
                   cg.quick(2 * l3p, 4, 2 * l3, 0, 0, 0);
  if (part2_0 != 0.0) {
    // part2_0 = part2_0 * PFK_ope_cal2(final_state,init_state,0,2,l3p);
    int num_p = wave.nl_vec[n12p][n3p][0][l3p];
    int num = wave.nl_vec[n12][n3][0][l3];
    part2_0 = part2_0 * wave.PFKglbar2_vec[num_p][num][0][0];
  }

  // lambda1 = 1; lambda2 = 1;
  double part2_1 = 0.0;
  int l_min = std::max(std::abs(l3 - 1), std::abs(l3p - 1));
  int l_max = std::min(l3 + 1, l3p + 1);
  if (l_min <= l_max) {
    double part2_1t = 0.0;
    for (int lbar = l_min; lbar <= l_max; lbar++) {
      double part2_1tt = 0.0;
      part2_1tt = (2 * lbar + 1) * cg.quick(2 * lbar, 2, 2 * l3p, 0, 0, 0) *
                  cg.quick(2 * lbar, 2, 2 * l3, 0, 0, 0) *
                  sixJ.quick(2, 2, 4, 2 * l3, 2 * l3p, 2 * lbar);
      if (part2_1tt != 0) {
        // part2_1tt = part2_1tt *
        // PFK_ope_cal2(final_state,init_state,1,1,lbar);
        int num_p = wave.nl_vec[n12p][n3p][0][l3p];
        int num = wave.nl_vec[n12][n3][0][l3];
        part2_1tt = part2_1tt * wave.PFKglbar2_vec[num_p][num][1][lbar];
      }
      part2_1t += part2_1tt;
    }
    part2_1 = part2_1t * 5.477225575051661; // sqrt(30)
  }

  // lambda1 = 2; lambda2 = 0
  double part2_2 =
      phase_B * std::sqrt(2.0 * l3 + 1) * cg.quick(2 * l3, 4, 2 * l3p, 0, 0, 0);
  if (part2_2 != 0.0) {
    // part2_2 = part2_2 * PFK_ope_cal2(final_state,init_state,2,0,l3);
    int num_p = wave.nl_vec[n12p][n3p][0][l3p];
    int num = wave.nl_vec[n12][n3][0][l3];
    part2_2 = part2_2 * wave.PFKglbar2_vec[num_p][num][2][0];
  }

  // cout<<" term B :  part1 = "<<part1<<"  part2_0 = "<<part2_0<<"  part2_1 =
  // "<<part2_1<<"  part2_2 = "<<part2_2<<endl;
  double result = part1 * (part2_0 + part2_1 + part2_2);
  return result;
}

double Ope_term::PFK_ope_cal0(const State_Jacobi &final_state,
                              const State_Jacobi &init_state) {
  double result = 0.0;
  int n12 = init_state.n12;
  int n3 = init_state.n3;
  int l3 = init_state.l3;
  int n12p = final_state.n12;
  int n3p = final_state.n3;
  int l3p = final_state.l3;
  if (l3 != l3p)
    return result;

  double part4 = 0.0;
#pragma omp parallel for reduction(+ : part4)
  for (int k12p = 0; k12p < wave.k_mesh; k12p++) {
    double part3 = 0.0;
    for (int k12 = 0; k12 < wave.k_mesh; k12++) {
      double part2 = 0.0;
      for (int k3p = 0; k3p < wave.k_mesh; k3p++) {
        double part1 = 0.0;
        for (int k3 = 0; k3 < wave.k_mesh; k3++) {
          part1 += wave.Pwave_vecwx[0][n3][l3][k3] * wave.F_reg_vec[k12][k3] *
                   wave.glbar_vec[0][l3][k3p][k3];
        }
        part2 += part1 * wave.Pwave_vecwx[0][n3p][l3p][k3p] *
                 wave.F_reg_vec[k12p][k3p];
      }
      part3 += part2 * wave.Pwave_vecwx[0][n12][0][k12];
    }
    part4 += part3 * wave.Pwave_vecwx[0][n12p][0][k12p];
  }
  result = part4;
  return result;
}

double Ope_term::PFK_ope_cal2(const State_Jacobi &final_state,
                              const State_Jacobi &init_state, int lambda1,
                              int lambda2, int lbar) {
  double result = 0.0;
  int n12 = init_state.n12;
  int n3 = init_state.n3;
  int l3 = init_state.l3;
  int n12p = final_state.n12;
  int n3p = final_state.n3;
  int l3p = final_state.l3;

  double part4 = 0.0;
#pragma omp parallel for reduction(+ : part4)
  for (int k12p = 0; k12p < wave.k_mesh; k12p++) {
    double part3 = 0.0;
    for (int k12 = 0; k12 < wave.k_mesh; k12++) {
      double part2 = 0.0;
      for (int k3p = 0; k3p < wave.k_mesh; k3p++) {
        double part1 = 0.0;
        for (int k3 = 0; k3 < wave.k_mesh; k3++) {
          part1 += wave.Pwave_vecwx[lambda2][n3][l3][k3] *
                   wave.F_reg_vec[k12][k3] * wave.glbar_vec[1][lbar][k3p][k3];
        }
        part2 += part1 * wave.Pwave_vecwx[lambda1][n3p][l3p][k3p] *
                 wave.F_reg_vec[k12p][k3p];
      }
      part3 += part2 * wave.Pwave_vecwx[0][n12][0][k12];
    }
    part4 += part3 * wave.Pwave_vecwx[0][n12p][0][k12p];
  }
  result = part4;
  return result;
}

void Ope_term::build() {
  int ch_size = jaco.channel.size();
  this->vec.resize(ch_size);
  Pgs_Bar pgs(ch_size, 30);
  // pgs.init()
  cout << "Begin to build Ope term :" << endl;
  for (int i = 0; i < ch_size; i++) {
    pgs.PrintBar(i);
    int size = jaco.state[i].size();
    vec[i].resize(size);
    // if(size<1)continue;
    // std::cout << "channel : " << i << endl;
    for (int x = 0; x < size; x++) {
      vec[i][x].resize(size);
      for (int y = 0; y < size; y++) {
        vec[i][x][y]=0.0;
      }
      for (int y = 0; y <=x ; y++) {
        int flag =-1;
        if (jaco.state[i][x].F > 0 && jaco.state[i][y].I > 0)flag=1;
        if (jaco.state[i][y].F > 0 && jaco.state[i][x].I > 0)flag=1;
        if (flag < 0)continue;
       //cout<<"channel : "<<i<<"\t x = " << x << "\t y = " << y << endl;
       double val = this->cal(jaco.state[i][x], jaco.state[i][y]);
        vec[i][x][y] = val;
        vec[i][y][x] = val;
      }
    }
  }
  pgs.EndBar();
  cout << "End to build Ope term :" << endl;
}

void Ope_term::print(int ch_num) {
  int size = vec[ch_num].size();
  for (int i = 0; i < size; i++) {
    for (int j = 0; j < size; j++) {
      if (vec[ch_num][i][j] != 0.0)
        cout << i << "\t" << j << "\t" << vec[ch_num][i][j] << endl;
    }
  }
}
void Ope_term::print() {
  for (int ch_num = 0; ch_num < vec.size(); ch_num++) {
    int size = vec[ch_num].size();
    cout << "channel num : " << ch_num << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        if (std::abs(vec[ch_num][i][j]) > 0.0001)
          cout << i << "\t" << j << "\t" << vec[ch_num][i][j] << endl;
      }
    }
  }
}

void Ope_term::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  for (int ch_num = 0; ch_num < vec.size(); ch_num++) {
    int size = vec[ch_num].size();
    file << "channel num : " << ch_num << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        if (vec[ch_num][i][j] != 0.0)
          file << i << "\t" << j << "\t" << vec[ch_num][i][j] << endl;
      }
    }
  }
  file.close();
}
